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We study the spectral function, S(k,u>) for the spin-1, one dimensional antiferromagnetic chain 
using a time-dependent density matrix renormalizaton group (DMRG) numerical method. We 
develop methods for extrapolating the time dependent correlation functions to larger times in order 
to enhance the frequency resolution. The resulting spectral functions are impressively precise and 
accurate. Our results confirm many qualitative expectations from non-linear a model methods and 
test them quantitatively. The critical wave-vector at which the single particle excitation emerges 
from the 2-particle continuum is estimated to be 0.23-7T — 0.24-7T. 



I. INTRODUCTION 



The 5=1 Heisenberg antiferromagnetic chain, 



if = J 2J 5,- • S J- 



(1) 



has been a subject of intense theoretical and experimen- 
tal study since Haldane's observation that it has an ex- 
citation gap above a singlet ground state to a triplet 
excited state, quite unlike the S=l/2 case. Much has 
been learned about the time independent properties us- 
ing a combination of analytic and numerical methods, 
in particular the density matrix renormalization group 
(DMRG).— Our knowledge of dynamical properties, while 
substantial, leaves more room for improvement. 2 Re- 
cently, DMRG has been significantly extended to al- 
low direct calculation of time dependence (tDMRG)&^£ 
and time-dependent correlation functions in particular. 
In this paper we will utilize tDMRG to calculate high- 
resolution spectral functions for the 5=1 chain for a 
broad range of momenta. Some of the numerical tech- 
niques developed here were briefly described in Ref. 0- 
The spectral function we shall focus on, in the isotropic 



S{k,w)= J2 



-ikj 



j=-oo 



dt e l " t (0|5 a (i)5 Q (0)|0) (2) 



where S{k,uo) = S aa {k,uj), with a = x, y, or z, and 
the subscripts on 5 indicate sites. The tDMRG method 
calculates the space-time dependent expectation values 
appearing in Eq. @ directly, and then one performs the 
Fourier transforms (FTs) in Eq. ||2J| to obtain S(k,iv). 
The crucial practical issue for this approach, which we 
discuss in detail in this paper, is dealing with the finite 
range of times available from a tDMRG simulation. We 
find that two different methods for extrapolation in time 
to increase the range of times used in the FTs are both 
very useful for increasing the frequency resolution. 

Much insight into the model can be obtained from an 
approximate mapping, based on the large S limit, onto 



the non-linear er-model (NLctM) with Lagrangian den- 
sity: 

£=^m$f-v 2 dj)% (3) 

where the lattice spin operators are represented as: 

s 3 ^—$xd t $+(-iy s $, (4) 

vg 

<f> = 1, v is the spin velocity and g ~ 2/s is the cou- 
pling constant, (j) IS t ne antiferromagnetic order parame- 
ter and I = (l/vg)4> X do4> is the uniform magnetization 
density. This low energy representation is only valid for 
wave- vectors near and tt. This field theory is known 
to have a singlet ground state in one space dimension, 
with a gap, A, to a massive spin-triplet of excitations. 
There are no bound states or any other single particle 
excitations besides this triplet. This implies that a single 
particle ((^-function) peak should appear in S(k,u>) near 
k = 7r at energy A but near k w the lowest excitations 
are a 2-particle continuum starting at 2A. Interaction 
effects mix <j> with (cf>) 3 [but not (<f>) 2 by symmetry] so 
that the theory also predicts a 3-particle continuum at 
fcsjir starting at 3 A, a 4-particle continuum at k w 0, et 
cetera. Continuity between t«0 and k w tt then implies 
a qualitative sketch like Fig. [T]for the region of non-zero 
spectral weight. In particular, the single particle peak 
must emerge out of the 2-particle continuum at some 
critical wave- vector k c . Integrability of the NLtrM allows 
for the calculation of exact form factors and hence predic- 
tions for the detailed shape of the 2-particle contribution 
to S(k, lS) near t«0 and the 3-particle contribution near 

Since the field theory is based on 5 — > oo and is only 
valid for k very close to and tt it is unclear how well 
any of these predictions should describe the S=l case. 
Various comparisons of NLcrM predictions with numer- 
ical results have been performed before, but these have 
necessarily focused primarily on equal time correlators 
and thermodynamic quantities. Experiments on quasi- 1- 
dimensional antiferromagnets have clearly confirmed the 
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FIG. 1: General features of the spectral function. The four 
and higher magnon continua are not shown. 



Haldane gap but the 2-particle nature of the small k ex- 
citations and the existence of a gap at k — ir from the 
single particle excitation at A to the bottom of the 3- 
particle continuum at 3A have not been confirmed and 
have led to some questioning of the validity of these field 
theory predictions. 

One of the purposes of this paper is to compare our 
tDMRG calculations of S(k,oj) with the NLcrM predic- 
tions. In addition, we will determine k c and examine the 
behavior near this wave- vector. 

In Sec. II we review the time dependent DMRG algo- 
rithm and discuss methods for extrapolating the time- 
dependent correlation function to longer times before 
Fourier transforming. In Sec. Ill we discuss properties 
of the single magnon excitation. In Sec. IV we discuss S 
near k ss and compare with the 2-magnon NLcrM pre- 
dictions. In Sec. V we study S near k rs n and compare 
with the 1-magnon plus 3-magnons predictions. In Sec. 
VI we study intermediate k in the vicinity of k c where 
the single magnon peak first emerges. Sec. VII contains 
our conclusions. 



II. TIME DEPENDENT DMRG AND 
EXTRAPOLATION IN TIME 

The time dependent DMRG method^i 5 - has been de- 
scribed elsewhere; here we describe some of the practical 
issues related to obtaining dynamical spectra. First, we 
describe the specific tDMRG algorithms used; second, 
time step and truncation errors; and third, extrapolation 
in time and windowing for time Fourier transforms. 

The first step of a calculation is to use ordinary DMRG 
to find the ground state of a finite open system to high 
accuracy, with typical lengths being L — 200 — 400. To 
avoid S = 1/2 end states, we put real S = 1/2 spins on 
the first and last sites. Let i be one of the two center 
sites. After obtaining the ground state <fi, we apply the 
operator or S? to create a state which is a mixture 
of excited states, ip(t — 0). We target and time evolve 



tp, and also, in order to minimize time step errors in 
correlation functions, <f>. As the time evolution occurs, 
we measure the correlation function S(J — i,t) when the 
tDMRG step is at site j, accumulating the function for 
a wide range of separations and times^ 

For the time evolution we use only the Suzuki- Trotter 
decomposition methods, which are very efficient for chain 
systems with only nearest neighbor interactions. We use 
two different variations. First, we use the original method 
of White and Feiguin^ which has second order Trotter 
errors. An advantage of this method is that every DMRG 
step has a time evolution bond operator applied, whereas 
Trotter methods splitting the links into even and odd 
groups 5 apply a time evolution operator on only half the 
steps of a sweep, resulting in roughly twice as much to- 
tal truncation error per unit time evolved. To obtain a 
measurement for a specific time, we perform the mea- 
surements in a half sweep without any time evolution. 
Thus, including measurements, we evolve in the repeated 
6 half-sweep pattern: evolve left-to-right, evolve right-to- 
left, left-to-right measurement half-sweep (without time 
evolution), evolve right-to- left, evolve left-to-right, right- 
to-left measurement half-sweep. During each half-sweep 
time evolution step, we evolve a time step r, so that mea- 
surements are available with a time step of 2r. Typically 
we use t = 0.1. (We set J = 1.) 

The second method we use is a fourth order Trotter 
methodf££ which with r = 0.1 virtually eliminates time- 
step errors, at the expense of more sweeps for a given 
time, and consequently larger accumulated truncation er- 
ror. We show below that the second order Trotter errors 
from the first method with a time step of r = 0.1 primar- 
ily lead to modest frequency shifts, shifting the Haldane 
gap by about 1%, for example. (The rate of decay in time 
of the correlation functions seems not to be strongly ef- 
fected.) For high accuracy studies it is more efficient to 
use the fourth order method rather than simply reducing 
r. The decomposition we use is£ 

e (A+S)r+0(r 5 ) _ gA0J gSflr gA(l-fl) } e S(l-28)r 

xe^-^e^e^* (5) 

with 9 = 1/(2 - 2 1 / 3 ) 1.35. Here A would represent, 
say, the odd bonds, and B the even, and each of the 
seven terms is applied in a half-sweep. Adding to this a 
measurement half-sweep, a total of eight half-sweeps are 
needed to evolve by r, with measurements available with 
a time step of t. 

Generally we specify a desired truncation error for each 
step, and vary the number of states kept m to achieve this 
truncation error. However, we also constrain m to be no 
larger than a specified miimit (typically 1000-2000), and 
no smaller than a minimum m m j n (typically 100-150). 
The purpose of m mm is to reduce the truncation error to 
near zero, at little cost, in steps where the state has very 
small entanglement. (These small entanglement steps oc- 
cur outside the "light cone" of the initial disturbance at 
site i, and time 0.) We specify a mn m it to avoid memory 



3 



limitations and to avoid a handful of steps (near site i) 
taking a large fraction of the computer time. 

The total accumulated truncation error, etot, summed 
over all DMRG steps since the start of the time evolution, 
is readily available and useful. At each step, a small 
part of each wavefunction is discarded; the truncation 
error is the magnitude squared of these small parts. In 
Fig. [2] we show that e to t gives a rough estimate of the 
typical errors to be expected in measurements at that 
time. The plot compares etot with the errors in S(x, t) 
for x = j — i = and x = 1. The errors in S(x,t) were 
estimated by running with two different truncation error 
parameters e, e = 2 x 10~ 10 and e — 4 x 10~ 10 . Wc 
also see that errors of order 10~ 5 are computationally 
feasible for times t < 10, but that the errors steadily 
grow with time. In this case the errors grow roughly 
linearly with time because a target truncation error per 
step was specified, with the upper limit TOii m ;t playing a 
small role. 
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FIG. 2: Comparison of the errors in measurements of S(x,t) 
and the total accumulated truncation error etot. Two runs 
were made using the 4th order Trotter decomposition, both 
with r = 0.1, with e = 2 x 10" 10 and e = 4 x 10" 10 . The 
values of AS were the difference in the results between these 
two runs. The value for etot was for the larger e run. 



To clarify further how this works in Fig. [3] we compare 
two runs, one a second order Trotter method, the other 
the more accurate run of Fig. [5] Both runs took compa- 
rable amounts of computer time. We define ro max to be 
the maximum number of states kept over all the steps in 
a half-sweep. Clearly m max < mi; m j t . In each half-sweep, 
the largest values of m were for steps near the center of 
the system, where the spin operator was applied. After 
a moderate time, for these steps near the center m was 
limited by mi; m it. The relatively small number of these 
steps made the effect of mu m i t on e to t small. The calcu- 
lation time for a step is proportional to rn 3 , so limiting 
m is important for efficiency. 

The tDMRG yields directly the space-time dependent 
correlation function directly. The signal initiated by the 
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FIG. 3: Maximum number of states kept m max versus time 
(upper- left curves), and total discarded weight etot versus 
time (lower right curves) for the systems of Fig. [5] The 
simulation adjusted m at each step to try to achieve a total 
discarded weight of e for that step, subject to a maximum 
m of either 1300 or 1800 (flat portions). The figure shows 
that in this case etot depends more strongly on e than on the 
maximum allowed m. 



application of the spin operator in the center of the sys- 
tem spreads out with time. Because the system has a 
finite correlation length, the correlation function is only 
non negligible within a range of about \vt\ of the center, 
where v is the maximum spin velocity. We always keep 
the maximum time i max of the simulation small enough 
so that the signal has not reached the edges at the end 
of the simulation. Hence there are essentially no spatial 
finite size effects. One can spatially Fourier transform 
(FT) x — > k for any k; the available values of k are con- 
tinuous, not discrete. On the other hand, the feasible 
^max is strictly limited by the available computer time, 
and the correlation functions decay slowly or not at all 
in time. Truncation of the signal following by an FT 
t — > u> would result in severe "ringing" . The standard 
approach is to multiply the signal by a windowing func- 
tion, which typically resembles a Gaussian centered at 
t = but which vanishes exactly at ±t max . A drawback 
of this approach is that most of the data gets "thrown 
away" , and the frequency broadening of the spectrum is 
large. 

To avoid this over-broadening, we have developed an 
alternative approach based on linear prediction.— Before 
the time-frequency FT, we extrapolate the time signal to 
long times using linear prediction. We then apply a broad 
window which does not throw away a significant amount 
of the original data. Linear prediction extrapolates a 
discrete equally spaced time series {yt} as 



(6) 



The coefficients dj are determined by the known data 
points {yi} by requiring that their prediction for each 
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point yi, based on yi- n ■ ■ ■ Ui-i, vary as little as possible 
from the actual value j/j, using a least-squares criterion. 
One finds that the dj are determined from correlation 
functions (yiUi+j), where the average is over i, and the 
principle computational work in determining dj is the in- 
verse ofanxn matrix. In our work we have used n = 20, 
so that the numerical work involved in the extrapolation 
is negligible. 
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FIG. 4: DMRG results for the real part of the onsite spin- 
spin correlation function S zz (x = 0,t) for the more accurate 
run of Fig. [2] The red curve represents DMRG data out to 
t = 22; the black curve is the DMRG data for t < 10 and a 
linear extrapolation for t > 10. 



In Fig. [J] we show the effectiveness of this extrapo- 
lation. For this run which had i max = 22, the data for 
t < 10 was extrapolated to longer times and compared 
with the DMRG results. The error in the extrapolation 
starts out small and grows reasonable slowly. Provided 
one does not rely on too long an extrapolation to try to 
achieve higher frequency resolution than the data sup- 
ports, this method performs much better than ordinary 
windowing methods. 

In Fig. [5] we compare the density of states, defined 
as N(ui) = S zz (x = 0,cl>), for the two systems of Fig. 
[3J After the linear extrapolation, the results were mul- 
tiplied by a Gaussian window exp[— t 2 / (20t 2 lax )] before 
Fourier transforming. On this scale the effects of the 
Trotter errors in the 2nd order data are almost not visi- 
ble. The results are considerably sharper than with the 
alternative method which did not use extrapolation. The 
sharp peaks are broadened square root singularities from 
the top and the bottom of the single-magnon dispersion. 
Above the top of the band, two- and three-magnon con- 
tributions to the spectrum are visible in a small tail. 
Using a larger window width would increase the reso- 
lution at the cost of increasing the likelihood of artifacts 
from the extrapolation. We will consistently use this win- 
dow width. Assuming the extrapolation is accurate, this 
means that our spectra should look like the exact spectra, 



but broadened by convolving with a gaussian 

CX P (-^LU 2 /W 2 ) 



(7) 



where the frequency resolution is W — (t max v 10)~ . 

Another approach is to fit the moderate time data to 
the correct asymptotic long-time form, in this case stem- 
ming from the square root singularities, to extrapolate 
to long times very accurately. The linear extrapolation 
asymptotically gives exponential decays in time, an in- 
correct assumption in this case, so for very long times 
the fitting method can be more accurate. A disadvan- 
tage of the fitting method is that it assumes one has 
some understanding of the results analytically. Another 
disadvantage is that the fitting process can take much 
more computer time than the linear prediction method, 
although still much less than tDMRG simulation itself. 
In cases where one does not know what sort of spectra to 
expect, one can first fit with the linear prediction method, 
and then guess an asymptotic form for fitting. 
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FIG. 5: Results for the density of states N(u>) for the two 
systems of Fig. [3] with linear extrapolation utilizing only the 
data out to t max = 20. The two curves are almost identical. 
The third curve used the 2nd order data but did not utilize 
extrapolation; instead, the data out to t ma x was multiplied by 
a simple finite window based on the cosine function. 



For the case of N(u>), we have leading singularities of 
the form 

a0{u - A)(w - A)~ 1/2 + W(n - - uj)- 1/2 (8) 

where f2 is the maximum in the single magnon dispersion 
relation near 7r/2 of about 2.725. This Fourier transforms 
to long time tails with leading terms of the form 

Aexp(-iAt)r 1/2 + Bexp(-iQt)t- 1/2 . (9) 

A more convenient form for fitting comes from the Fourier 
transform integral identities 



due' 1 " 1 8(uj - b)e- a ^-V {ui - b) 9 

> 

= r(l+ g)e- M {a + it)- 1 -3 (10) 



5 



and 

/DO 
dtue- lut 0(b - uj)er< b -^ (b-U))9 
-OO 

= r(l + g)e- M {a-ity 1 -a. (11) 

These identities are useful because the frequency expres- 
sions have only a single one-sided singularity, and con- 
venient well-behaved Fourier transforms. Eq. (flTj]) is 
useful for describing a singularity on a lower edge, while 
Eq. (fTT|) describes an upper edge. Thus to fit to the 
tDMRG data for the Fourier transform of N{u>) : we use 
the asymptotic form 

Aexp(-iAt)(a + ity 1/2 + B exp(-int)(b-ity 1/2 . (12) 

Fitting the accurate 4th order data over the time range 
10 — 20 with this form, we find that the fit within this 
range matches the data very accurately, with a typical 
absolute deviation of about 2 x 10 -4 , or a relative error 
of about 10~ 3 . 

Using the fitting parameters and asymptotic form, one 
extends the data to large times. We did not use the 
identities Eqs. (fTU)) and (fTl"j) to help perform the Fourier 
transform (although this might be convenient); we simply 
used a fast Fourier transform over a very large range of 
times (e.g. -10000 < t < 10000). Over the fitting region 
in time a smooth transition is made from the data for 
small times to the fit for large times. The results of this 
procedure for N(uj) are shown in Fig. [5] The fits allowed 
both A and f2 to vary; the result for the gap for the time 
range 10 - 20 was A = 0.4104327, accurate to 4 digits 
(see next section) . The resulting spectra, based on fitting 
over different time ranges, appears to be accurate to the 
line width in the figure. 

III. SINGLE MAGNON RESULTS 

Fourier transforming the DMRG S(x,t) data x — ► k, 
we obtain S{k, t) for any k. The time Fourier trans- 
form gives spectra with (potentially) single magnon and 
multimagnon contributions. As we discuss below, for 
fc c < fc < 7r, with k c w 0.237T — 0.247T, a single magnon 
delta function peak is present, plus multimagnon con- 
tinua at higher frequencies. In this section we focus on 
this well-defined single-magnon mode. 

First consider the band minimum, at k = tt, where the 
excitation energy is the Haldane gap. The most accu- 
rate method to determine the Haldane gap is still ground 
state DMRG, where we have used the same method as 
in Ref. US but with up to m = 500 and L = 400, 
to determine the Haldane gap to very high accuracy, 
A = 0.41047925(4). (The end coupling used to make 
the lowest excitation have k = ±tt was J cn( j = 0.50865, 
compared to 0.5088 in Ref. 0.) 

In Fig. \7\ we show results for S(tt, ui) near to = A. Lin- 
ear extrapolation plus Fourier transforming as described 
above give rather narrow Gaussian-shaped peaks. The 
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FIG. 6: Results for the density of states N(u>) from fitting to 
the asymptotic form Eq. (|12[) over the time range indicated 
in the caption, and then extending to very large times with 
the fit before Fourier transforming. The two curves overlap, 
indicating that the error in the determination of N(uj) is less 
than the line width. 



2nd order peak is narrower because of a larger i max , and 
it is shifted from the exact result because of Trotter er- 
ror. Because two and three magnon contributions are 
very weak and separated in frequency from A, a least 
squares fit to a pure exponential is almost identical to 
the maximum of the broadened peaks. In fact, neglect- 
ing any Trotter shifts, either the maximum or the fit 
frequency provide much more accurate determinations of 
the exact magnon energy than the peak widths would 
indicate. For a set of fc's spaced 0.017T apart we have fit 
the time data either to a pure complex exponential (for 
k > 1) and for k c < k < 1 to a complex exponential plus 
an asymptotic form describing the near-by two-magnon 
edge. These latter more complicated fits are discussed in 
Section VI. The frequency of the exponential term deter- 
mines the dispersion e(k) for k > k c . With this fitting 
approach, a larger i max is not very important compared 
to the Trotter error, so we utilize the 4th order data. The 
corresponding result for e(7r) is 0.41050; the error is only 
a few times 10~ 5 . For smaller k the error is expected to 
be larger because the multimagnon continuum is larger 
relative to the single magnon peak, and the continuum is 
closer in frequency, but the errors for k > k c are probably 
no bigger than 10~ 3 . 

In Fig. [5] we show e(k) from the 4th order run, along 
with an analytic fit motivated by the NLer model. The 
massive triplet excitations of the NLcrM have the rela- 
tivistic dispersion relation: 

e (fc) = \/A 2 +v 2 k 2 . (13) 
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FIG. 9: Comparison of the single magnon dispersion from 
DMRG and the relativistic approximation, Eq. 1131 



FIG. 7: Results for S(tt, u) for the two systems of Fig. [3] The 
Gaussian shaped curves come from linear extrapolation and 
Fourier transforming. The vertical lines represent delta func- 
tions coming from a least squares fit of Ae 1 ^ 1 to S(7r,t). The 
4th order fit frequency and the "exact" ground state DMRG 
result are indistinguishable in this plot. 




FIG. 8: Single magnon dispersion. The curve labeled 
"DMRG" comes from the fit of a pure exponential to the 
DMRG S(k, t) data. The curve labeled "Fit" is the analytic 
expression Eq. (|14|) . The final curve is the two magnon band 
minimum at 2e(7r — k/2). The inset shows the region near 
fc c , where the magnon line enters the two magnon continuum. 
The first two curves are not very meaningful well below fc c , 
since the there is no single magnon delta function. 



Here the momentum k of the NLcM is k — it for the spin 
chain. In the NLctM, k can take any real value. Of course, 
in the spin chain, crystal momenta lie in the Brillouin 
zone, |fc| < 7r. This discrepancy limits the validity of 
the NLerM, especially when we consider multi-particle 
excitations, eo(fc) is only an approximation to the exact 



single magnon dispersion relation, e(k). We expect a 
perfectly stable single magnon excitation to exist for k < 
7T— k c with this dispersion relation. Strictly speaking s(k) 
is not defined for k > ir — k c . The "Fit" curve shown in 
Fig. [HI is based on the expression 



e(fc) « A 



\ 



1 + a n {l — cos[?i(7r — k)]}. 



(14) 



with the parameters a n given in Table 1. The gap, at 
k — 7r, from the data used in the fit is A = 0.410504; one 
could also use the more accurate value A = 0.41047925. 
e(k) goes through a maximum of 2.72551 at k w .47677 
and has the value 1.96 at k « 0.237T near where the single 
magnon excitation becomes unstable. There is an inflec- 
tion point (d 2 e/dk 2 — 0) at h ln » 0.8687T. As shown 
in Fig. [S] e(fc), agrees quite well with the Lorentz in- 
variant approximation, eo(fc), for k < O.ln , and rea- 
sonably well for k < 0.2ir, with v w 2.472. Also shown 
in Fig. [5] is the two magnon band minimum, which for 
\k\ < 2(tt - k in ) fn 0.265tt is given by 2e(7r - k/2). Near 
k 0, this is approximately 



E (fe)-> A + 



{vkf 
2A 



(15) 



The intersection of the two magnon minimum and the 
single magnon dispersion line determines k c , as show in 
the inset of Fig. [51 

The single magnon dispersion can be used to construct 
the multimagnon band minima and maxima, assuming 
that the magnon-magnon interactions are negligible, and 
adding the single magnon energies. The resulting bands 
are shown in Fig. 1101 In this construction, it was as- 
sumed that the single magnon line stops abruptly at 
0.237T ps k c . This resulted in slope discontinuities vis- 
ible in the two magnon band. However, at k c , as we 
discuss below, the sharp " single magnon peak" can be de- 
scribed as either a single magnon or two magnon feature. 
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FIG. 10: Single magnon line (black solid line), two magnon 
band (horizontal line fill, red) , and three magnon (dotted fill, 
green) bands. 



This ambiguity effectively blurs the distinction between 
the two and three magnon bands near the slope discon- 
tinuities. If one changed the construction to include the 
broadened "single magnon peak" to values of k below k c , 
the two magnon bandwidth would be broadened near the 
slope discontinuities. This suggests that, for example, at 
k = 7r one might expect to see the two magnon band 
maximum as a visible feature of the spectrum, but one 
might not see anything for the two magnon minimum. 
(The spectrum at k = 7r is shown in Section V.) 

The amplitude of the single magnon peak is approxi- 
mated well byii 



S(q,t = 0) 



vZ 



(16) 



with Z = 1.26. A related quantity is the fraction of the 
spectral weight in the multimagnon continua, given by 



/ = 



S(q,t = 0)-A 
S(q,t = 0) 



(17) 



where A is the amplitude of the single magnon peak. 
This quantity is shown in Fig. 111! It is interesting that 
/ shows non monotonic behavior with k, with minima 
near k = it/ 2 and k = tt. This non monotonic behav- 
ior is roughly correlated with the gap between the single 
magnon line and the lowest multi-magnon band mini- 



TABLE I: Coefficients in the fit of the single magnon disper- 
sion relation, Eq. (|14[) . 



1 1.96615 

2 21.20162 

3 -1.61279 

4 -0.04766 

5 0.02407 



mum. Naively, one might expect that a nearby contin- 
uum has an easier time than a faraway one does in taking 
spectral weight from the single magnon line. 




0.25 



0.5 0.75 
k/jc 



FIG. 11: Fraction of the spectral weight in the multimagnon 
continua for k > k c , shown by the solid black line. The gap 
between the single magnon line and the lowest multimagnon 
band is shown by the dashed red line. 



IV. k as 

We now discuss the properties of the spectral function 
for small k. That S oc k 2 as k — > can easily be proven to 
be exactly true using the fact that J^ - 5||0) =0 (singlet 
ground state) and Taylor expanding S(k, w) in Eq. © to 
second order in k. On general symmetry grounds, based 
on the NLtrM, we expect that near k as 0, S will con- 
tain only multi-particle continua corresponding to even 
numbers of bosons. At a finite small k, the lowest en- 
ergy 2-magnon state with total momentum k is one in 
which each magnon has momentum tt + k/2. Therefore, 
since the dispersion is even about n, the bottom of the 
2-magnon continuum should be exactly at 2e(7r — k/2). 
This remains true up to k — 2(tt — fc in ) as .2657T (recall 
k[ n as 0.868-7T is the inflection point). For a range of en- 
ergies, at small enough k, the only possible excitations 
have 2 magnons. This is true up to u> — 4e(fc/4), for 
k < (tt - k in ). 

A simplified, "mean field" version of the NLcM is a 
free massive boson model with Lagrangian: 



1 

2^ 



[(^) 2 - w 2 (a^) 2 -A 2 (0) 2 ], 



(18) 



and no constraint on <j). Expanding <f> in boson creation 
and annihilation operators one finds the free boson result 



S (k,u)) 



k 2 y/LU 2 - (kv) 2 -4A 2 

v[uj 2 — (kv) 2 } 3 / 2 



8{uj 2 - (vk) 2 



4A 2 



As expected from general principles, this vanishes 
quadratically as k — > 0, and also vanishes below the 2- 
magnon threshold, Wth = 2^ A 2 + (vk/2) 2 . The exact 
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2-magnon expression for S(k,u>), in the NLerM is known 
exactly and can be written: 



S 0a (k,Lu) = S {k,Lu) 



7T 4 1 + (6»/tt) 2 /tanb.0/2 



64 1 + (0/2tt) 2 V 6/2 
where the rapidity, 0, is defined by: 

6 = 2cosh- 1 [(w 2 - v 2 k 2 )/(4A 2 )]. 



(20) 
(21) 



2nd order 

4th order 

-3/2 




40 



50 



FIG. 12: Comparison of \S(k, t) | for k = n/W for two different 
runs, with 2nd order and 4th order Trotter decomposition 
runs of Fig. 2, and the asymptotic form expected analytically. 



Just above wth, both Eqs. (|19| and (|20|) rise as 
{lo — Wth) 1 / 2 , which would lead to the asymptotic time 
behavior 



S(k,t) 



-iuit^t ±— 3/2 



(22) 



In Fig. [12] we compare \S(k,t)\ for a typical small value 
of k with At~ 3 / 2 , finding good agreement for large t with 
the empirical parameter A = 0.6. Using the linear predic- 
tion method, we Fourier transformed the S(k,t) DMRG 
results for both runs, and in Fig. Q2] they are compared 
to Eqs. (flU]) and (f2T)]) . They look qualitatively similar. 
In particular, the threshold singularity appears the same 
and they have peaks at similar frequencies. However, the 
peak is about a factor of 2 larger in the DMRG data than 
in the NLctM. Furthermore, the field theory results drop 
off much more slowly at large w. This latter feature is 
to be expected since magnons of arbitrarily high momen- 
tum are included in the field theory while there is a cut 
off at the Brillouin zone boundary in reality. 

A simple way to improve the high frequency behavior 
in the free boson approximation is to replace the rela- 
tivistic model of Eq. (fP5|) , by a Hamiltonian: 



H 



—k + e- 2 (k)<j)k ■ <t>-k], 



(23) 



Free boson model 
Nonlinear a model 
4th order 
2nd order 




FIG. 13: S(k,u>) for k = 7r/10 for two different runs, with 
different accuracies and different total times, both Fourier 
transformed using linear extrapolation. For comparison, two 
analytic results, based on Eqs. (|19|l and (I20|l . are shown. 



numerically determined single magnon dispersion rela- 
tion. The appropriate form of the small k spin operators: 



S k 



k' 



x n 



k—k' 



(24) 



is determined by the requirement that the spin commu- 
tation relations are obeyed and that Sq commute with 
the Hamiltonian. This changes So{k,Lu) to: 



Sq 



[e(fc') - e(k - k')} 2 



2e(k')e(k - k')\e'{k') - e'(k - fc')l ' 



(25) 



where e(k) is the exact (numerically determined) disper- 
sion relation and e'(k) denotes its derivative, k' , the mo- 
mentum of one of the 2 magnons, is determined from u> 
and k by energy-momentum conservation: 



e(fc') + e(k - k') =U3. 



(26) 



For low enough energy, there is only one pair of solutions 
to Eq. (|2"6"]) . and one element of the pair should be chosen 
in evaluating Eq. (|25[) . The resulting improvement of the 
free boson result is shown in Fig. [14] As expected, there 
is little change near the peak and threshold, but the high 
energy tail is cut off. (Here we restricted each boson to 
have |fc| < 7r/2.) 



In addition to the single magnon mode, which has most 
of the spectral weight near k = tt there is also a contri- 
bution from 3, 5, . . . magnons. Of these, the largest is 
expected to be the 3-magnon contribution. The lower 
threshold for this, corresponding to each magnon hav- 



where n£ is canonically conjugate to ^ and e(k) is the ing momentum fc/3 is at 3e(fc/3) ps 3W A 2 + (vk/3) 
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FIG. 14: S(k,u>) for k = n/10, comparing the near-exact 
DMRG results with the modified free boson result. 



At small k the exact lower threshold approaches 3A + 
(vk) 2 /(6A). Note that the presence of a multi-magnon 
contribution is a consequence of inter-magnon interac- 
tions; it vanishes for the non-interacting model of Eq. 
(fTHj) . The exact 3-magnon form factor is known for the 
NLcrM and the resulting 3-magnon contribution to S can 
be expressed in terms of an elementary integral. The re- 
sult is compared to our DMRG results in Fig. [T5] Again 
there is a qualitatively similarity with a peak at a simi- 
lar energy but now the NLcrM peak is about 3 times too 
low and there is far too much spectral weight at high en- 
ergies. The total spectral weight in the 3-magnon peak 
compared to that in the single magnon is found from 
DMRG to be 2.7%. The lower edge of the multimagnon 
band is given by the three magnon edge, c.f. Fig. 1101 
The band has a sharp dropoff at the two magnon band 
maximum. The three magnon band above that is rather 
small. As discussed earlier, one does not expect a sharp 
feature for the two magnon minimum (nominally near 
w/A « 9), and none is visible. 



VI. k « k c 



FIG. 15: S(k, ui) for k — n in the multimagnon frequency 
regime for two different runs. Also shown are the results from 
the NLcrM model. The vertical line near uj/A ~ 13 is two 
times the maximum of the single magnon dispersion, roughly 
locating the top of the two magnon band. The tiny bump 
near uj/ A = 2 is an artifact of the FT of the much larger 
single magnon peak at ui/A — 1, while we believe the small 
tail above ui/A ~ 13 is real. 




FIG. 16: S(k, uj) near k c » 0.23tt - 0.24tt for the 4th order 
run using the linear prediction method. The vertical line as- 
sociated with each k is the two magnon lower band edge. 



A remarkable feature of S(k,Lu) which is completely 
missed by the NLcrM approach is the merging of the 
single particle peak into the 2-particle continuum at 
k = k c fa 0.23?r - 0.24tt. Results for S(k,u) from the 
linear prediction method near k c are shown in Fig. 1161 
Above k c , one sees the separate single magnon peak, 
broadened by the finite run time and Fourier transform. 
By k = 0.27T the peak has disappeared, and one sees a 
characteristic small k line shape. Close to k c the broad- 
ening from the finite maximum time obscures the details 
of the spectrum. 

The fitting method is very useful to capture the behav- 
ior of S(k,u>) near k c more accurately. Just above k c we 
expect a combination of a single magnon S function peak 



and a two-magnon continuum similar to that at small k. 
To fit the time data as accurately as possible, we assume 
that the two magnon band starts exactly at the expected 
threshold w th = mim./ e(ir-k/2+k')+e(ii~k/2-k'). The 
relevant values of e(k) used to determine u>th are far from 
k c and a simple exponential fit determines the magnon 
peak location very accurately. We also assume that the 
small k threshold behavior (uj — wth) 1 ^ 2 applies. Includ- 
ing, in addition, the next expansion term (uj — w t h) 3 ^ 2 , 
we utilize the fitting form (cf Eq. (|12p) 

Aexp(-iujt) + Bexp(-iuj th t)(b + it)^/ 2 

+Cexp(-iuj th t)(c + it)- 5 / 2 . (27) 
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Results from this fitting followed by Fourier transforming 
are shown in Fig. Q~7] The fits deviated from the data 
over the range t = 10 — 20 typically by a few times 10~ 5 , 
and the magnitude of the data points fitted to was typi- 
cally near 0.1 — an excellent fit, making a convincing case 
that the assumed asymptotic form is correct and that the 
results for S(k,u>) are very accurate. The peak locations 
from these fittings were used in the determination of the 
dispersion relation of Section III near k c . 
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FIG. 17: S{k, u) near k c w 0.23tt - 0.24tt for the 4th order 
run using fitting to the asymptotic time decay. 



One can view the split-off of the single magnon peak 
from the two particle continuum at k c in two different 
ways. First, in an approach motivated by the small k 
two-magnon nature of the spectrum, one can regard the 
splitoff as due to the formation of a sharp 2-magnon 
bound state for k > k c ^ Second, if we imagine that 
k c is "large" so that the two magnon description is inap- 
propriate, then we can regard the single magnon peak as 
surviving for k near but less than k c , but with a broad- 
ening caused by decays into two magnons. Here we will 
consider in more detail these two pictures. 

Consider first the two-magnon bound state picture. 
Since the quantum numbers of the single magnon state 
with S z = 1 and a two magnon bound state composed 
of a S z = and a S z = 1 magnon are identical, we 
are free to regard the excitation which splits off in either 
way. Well above k c the two-magnon bound state pic- 
ture is clearly not very useful, since its formation would 
imply large magnon-magnon interactions which are not 
otherwise observed. The NLcrM does not have any bound 
states. Furthermore, in any Lorentz invariant theory, in- 
creasing the centre of mass momentum can never lead 
to bound state formation. As was discussed in Ref. [ToL 
since the bosonic magnon excitations form a triplet, and 
any excitation produced from the singlet ground state by 
the spin operators must also be a triplet, it follows that 
the bound state wave-function must be antisymmetric in 
its spatial coordinates. In Ref. [l(| it was found that the 
effective magnon-magnon interaction is attractive in the 



antisymmetric channel. On the other hand, it is appar- 
ently repulsive in the symmetric channel. This is related 
to the BEC picture of the transition at a critical magnetic 
field 1 ^ where the Zeeman energy equals A. 

Even in one dimension, an arbitrarily weak attraction 
does not produce a bound state in the antisymmetric 
channel (although it does produce one in the symmetric 
channel.) It was observed in Ref. [lOj that as k is increased 
from zero, the momenta of the two magnons forming the 
bound state near the threshold u> m i n (k), which are near 
7r + k/2, approach the inflection points ki n ~ 1.1317T. To 
study the bound state near the threshold, we can expand 
the dispersion relation, e, near tt + k/2. The momenta 
of the two bosons are: tt + k/2 ± q. Expanding the total 
kinetic energy in powers of q gives the effective kinetic 
energy : 

P (4) 

T^2, + ^q 2 + — 4 + ..., (28) 

where e and its derivatives are evaluated at 7r + k/2. We 
see that the effective mass for the centre of mass motion 
is given by: 

±- = e(V(n + k/2). (29) 

(1/to) vanishes at the inflection point, k — > 2km ~ 2ir sa 
0.262-7T and the effective kinetic energy becomes quartic. 
At ki n , the coefficient of the quartic term is e( 4 )/12 ~ 
3.341 > 0. For a quartic kinetic energy, an arbitrarily 
weak attraction leads to a bound state, in both symmet- 
ric and antisymmetric channel. (For the antisymmetric 
case this can be seen by considering a trial wave-function 
of oc xe~~ x /( 2w >. The potential energy can be assumed 
to be everywhere less than a square well of depth vq 
and width a. For w ^> a, the potential energy is less 
than a quantity oc —vo(a/w) 3 , while the kinetic energy is 
oc (l/u>) 4 . For small vq, this has a negative minimum at 
w oc 1/vq 3> a, proving the existence of an antisymmetric 
bound state with binding energy oc v 4 .) This argument 
implies that k c < 2\n — ki n \ k, .262tt, since for any attrac- 
tive potential, as the effective mass diverges the bound 
state will eventually form. It is interesting to note that 
our DMRG estimate of k c ~ 0.23-7T — 0.24-7T is only very 
slightly less than 2|7r — fc, n |, suggesting that the attrac- 
tive interaction between magnons is weak. In addition, 
the splitting of the peak from the two magnon continuum 
varies as (to — mo) 2 (where mo is the mass at k = k c ) 
and thus as [k — k c ) 2 within this bound state picture. 

Now consider the second picture of the split-off, that 
of a magnon entering the continuum, but surviving in a 
broadened form near k c . In this case one would expect 
for k > k c the single magnon peak and the two magnon 
continuum would vary independently with k and that the 
splitting of the peak from the two magnon continuum 
would be linear in k — k c . In Fig. 1181 we show this 
splitting near k c . Indeed, the splitting appears to be 
linear in k — k c . Below k Cl the broadening grows very 
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rapidly. If the splitting is quadratic in k — k c , it must be 
so only very close to k c . 
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FIG. 18: Splitting between the single magnon peak and the 
bottom of the two magnon band near k c determined from a fit 
of the time dependent DMRG data. For uj > 0.24, the data 
was fit well assuming the single magnon peak was a delta 
function; below a better fit was obtained assuming the peaks 
was broadened as a Gaussian. The "error bars" indicate the 
width of the Gaussian. 



The behavior of the splitting versus k — k c favors the 
second picture of the split-off. However, we cannot rule 
out that the two-magnon bound state picture applies very 
close to k c . Note that within the two magnon bound state 
picture, there does not seem to be a compelling reason 
for more than a small fraction of the spectral weight to 
appear in the bound state. In fact, as shown in Fig. Ill) 
about 50% of the spectral weight appears in the single 
magnon peak near k sw 0.25, but the weight in the peak 
is rapidly falling as k is decreased. 

It is also interesting to examine the line shape for k 
slightly less than k c . In order to compare line shapes 
for different k's, in Fig. [19] we have shifted the curves 
to make the two-magnon thresholds identical, and have 
scaled them to make them identical at the arbitrary point 
uith + 0.25. These curves were made using the linear pre- 
diction method. The curve for k — 0.2ir shows a sharp 
resonance persists below but near k c . This resonance 
disappears by the time k is reduced to 0.14-7T. 



has proved to be an extremely effective method for cal- 
culating spectral functions for the 5 = 1 chain. We have 
been able to study fine details of the spectra with much 
greater resolution and accuracy than with any previous 
method. In comparing with free boson and nonlinear 
sigma model predictions for features of the spectra near 
k = and k — n, we find good qualitative agreement, 
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FIG. 19: S(k, oj) near and below fc c « 0.23vr - 0.24tt for 
the 2nd order run. Each curve is shifted by the two-magnon 
threshold energy 2e(fc/2), and scaled by an arbitrary factor 
to make S(k,u — Wth = 0.25) identical in each of the three 
curves. 



but quantitative disagreements in the overall magnitude 
of the spectrum and in the high frequency tails. Our re- 
sults near k c where the single magnon peak enters the 
two magnon continuum are better described in terms of 
a single magnon exhibiting decay and scattering below 
k c rather than viewing the single magnon peak as the 
formation of a two magnon bound state above k c . 



VII. CONCLUSIONS 

The combination of time dependent DMRG and ex- 
trapolation of the time dependent correlation functions 
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